%plots perturbation population response
%plots population resonse on running onset 

ii=figure;
oo=figure;
mean_act=[];
mean_act_run=[];
sem_act=[];
for site=[1:3,5:size(proj_meta,2)]%1:size(proj_meta,2)
no_time_points = size(proj_meta(site).rd,2);

act=[];
type_sessions={'f','p','d'};
for tp=1:no_time_points
    act=[];
    for lyr=1:size(proj_meta(site).rd,1)
        act=[act; proj_meta(site).rd(lyr,tp).act];
    end
    %     act=act(max(act')>2,:);
    sess_start=proj_meta(site).rd(lyr,tp).nbr_frames;
    sess_start=cumsum(sess_start);
    sess_start=[1 sess_start+1];
    velM = proj_meta(site).rd(lyr,tp).velM_smoothed;
    velP = proj_meta(site).rd(lyr,tp).velP_smoothed;
    
    ac_sess=find(strcmp(proj_meta(site).rd(lyr,tp).session,type_sessions(1))==1);
    indices=[];
    for kk=1:length(ac_sess)
        indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
    end
    
    velM_r = velM>0.001;
    start_perturb=strfind((proj_meta(site).rd(lyr,tp).ps_id>1),[0,1]);
    start_perturb=intersect(indices,start_perturb);
    ind_perturb=[];
    for kkk=1:length(start_perturb)
        if start_perturb(kkk) > 5
            if sum(velM_r(start_perturb(kkk)-10:start_perturb(kkk)+50))>52
                ind_perturb = [ind_perturb start_perturb(kkk)-10:start_perturb(kkk)+50];
            end
        end
    end
    temp2=act(:,ind_perturb);
    temp2= reshape(temp2',61,numel(temp2)/61)';
    mean_act(site,1:61)=mean(temp2,1);
    
    %% running onset
    % first fb session - filter out the running and stim moving
    ac_sess=find(strcmp(proj_meta(site).rd(lyr,tp).session,type_sessions(2))==1);
    indices=[];
    ind_run=[];
    for kk=1:length(ac_sess)
        indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
    end
    velP_r = find(velP<0.001==1);
    temp_start_run = velM_r(1:end-2)+velM_r(2:end-1)+velM_r(3:end);
    temp_start_run(temp_start_run>1.5)=3;
    temp_start_run(temp_start_run<=1.5)=0;
    temp_start_run = strfind(logical(temp_start_run),[0,1]);
    start_run=intersect(indices,temp_start_run);
    start_run=intersect(start_run,find(velP_r==1));

    for kkk=1:length(start_run)
        if start_run(kkk)>10 && start_run(kkk)+30<sess_start(end)-1
        if sum(velM_r(start_run(kkk):start_run(kkk)+30))>28 && ...
            sum(velP_r(start_run(kkk):start_run(kkk)+30))<28
                ind_run=[ind_run start_run(kkk)-10:start_run(kkk)+50];
        end
        end
    end
    % for dark
    ac_sess=find(strcmp(proj_meta(site).rd(lyr,tp).session,type_sessions(3))==1);
    indices=[];
    for kk=1:length(ac_sess)
        indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
    end
    start_run=intersect(temp_start_run,indices);

    for kkk=1:length(start_run)
        if start_run(kkk)>10 && start_run(kkk)+30<sess_start(end)-1
        if sum(velM_r(start_run(kkk):start_run(kkk)+30))>28
            ind_run=[ind_run start_run(kkk)-10:start_run(kkk)+50];
        end
        end
    end
    temp2=act(:,ind_run);
    temp2= reshape(temp2',61,numel(temp2)/61)';
    mean_act_run(site,1:61)=mean(temp2,1);
    
end

mean_act(site,:)=bsxfun(@rdivide,mean_act(site,:),mean(mean_act(site,1:5),2));
mean_act_run(site,:)=bsxfun(@rdivide,mean_act_run(site,:),mean(mean_act_run(site,1:5),2));

figure(ii);
colors=[0 0 1;0 0 1;1 0 0;1 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;];
set(0,'DefaultAxesColorOrder',colors(1:no_time_points,:));
subplot(4,2,site);
plot(mean_act(site,:)')
% plot(mean_act');
title([proj_meta(site).animal ' - ' num2str(proj_meta(site).ExpGroup(1))]);

figure(oo);
colors=[0 0 1;0 0 1;1 0 0;1 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;];
set(0,'DefaultAxesColorOrder',colors(1:no_time_points,:));
subplot(4,2,site);
plot(mean_act_run(site,:)')
% plot(mean_act');
title([proj_meta(site).animal ' - ' num2str(proj_meta(site).ExpGroup(1))]);

end
figure;
hold on,
mean_act=mean_act([1:3,5:7],:);
mean_act_run=mean_act_run([1:3,5:7],:);
% plot(mean_act','r')
% plot(mean_act_run','b')
% plot(mean(mean_act,1),'r.')
% plot(mean(mean_act_run,1),'b.')
% plot((mean(mean_act,1)+std(mean_act)/sqrt(6))','k+','linewidth',2)
% plot((mean(mean_act,1)-std(mean_act)/sqrt(6))','k+','linewidth',2)
% plot((mean(mean_act_run,1)-std(mean_act_run)/sqrt(6))','k+','linewidth',2)
% plot((mean(mean_act_run,1)+std(mean_act_run)/sqrt(6))','k+','linewidth',2)

aa=figure;hold on
ac=(mean(mean_act)-1)*100;
ac_sem=((std(mean_act)/sqrt(6)))*100;
ac2=(mean(mean_act_run)-1)*100;
ac_sem2=(std(mean_act_run)/sqrt(6))*100;
area([-10/10:1/10:50/10 50/10:-1/10:-10/10],[ac + ac_sem fliplr(ac-ac_sem)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');
area([-10/10:1/10:50/10 50/10:-1/10:-10/10],[ac2 + ac_sem2 fliplr(ac2-ac_sem2)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');

plot(-10/10:1/10:50/10,ac,'Color',[244/256 164/256 96/256],'Linewidth',3)

plot(-10/10:1/10:50/10,ac2,'Color',[0 0 0],'Linewidth',3)
xlim([-1.1 5.2])
ylim([-2 7])
xlabel('time [s]')
ylabel('\DeltaF/F [%]')

